Precise Simulation of Near-critical Fluid Coexistence 
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We present a novel method to derive liquid-gas coexisting densities, p^(T), from grand canonical 
simulations (without knowledge of Tc or criticality class). The minima of Ql = / {'m'*') l in an 

LxLxL box with m — p — {p)l are used to generate recursively an unbiased universal finite-size 
scaling function. Monte Carlo data for a hard-core square-well fluid and for the restricted primitive 
model electrolyte yield to ±1-2% of pc down to 1 part in 10^-10'^ of Tc (and confirm well Ising 
character). Pressure mixing in the scaling fields is unequivocally revealed and indicates Yang- Yang 
ratios = —0.044 and 0.26 for the two models, respectively. 

PACS numbers: 64.60.Fr, 02.70.Rr, 05.70.Jk, 64. 70. Fx 



Determining phase boundaries, critical points, and uni- 
versality classes for various models that lack a clear sym- 
metry has presented a serious difficulty in computer sim- 
ulations [1,2]. To tackle this problem, understanding 
scaling behavior in systems of finite size is crucial. How- 
ever, as recently stressed [3], an important issue arises for 
asymmetric fluid criticality, even in the thermodynamic 
limit, namely, the potential presence of a Yang- Yang 
anomaly, in which the second derivative of the chemi- 
cal potential, fia-iT), on the gas-liquid phase boundary 
diverges when the critical point, Tc, is approached from 
below. To describe a Yang- Yang anomaly requires pres- 
sure mixing in the scaling fields [3-5]. This also gen- 
erates a term varying as {t]"^^ [with t={T ~ Tc)/Tc] in 
the gas-liquid coexistence diameter, that dominates the 
previously recognized |t]^~" term [6] and further distorts 
coexistence curves near criticality. 

Our aim here is to show how coexistence curves may be 
estimated precisely and reliably near asymmetric critical 
points using grand canonical simulations, and to check 
our current understanding of scaling in such cases [4,5]. 
It transpires that a finite-size scaling analysis at Tc also 
elucidates pressure mixing and allows us to measure its 
strength using simulation data. 

Figure 1 presents our estimates of Apoo {T) = p^ — p^' , 
the density discontinuity across the phase boundary, for 
a hard-core square-well (HCSW) fiuid and for the re- 
stricted primitive model (RPM) electrolyte, where p^(T) 
and p~ (T) are the coexisting densities of liquid and va- 
por. The crosses represent new estimates obtained as 
explained below, while the open circles were derived pre- 
viously directly from the observed double-peaked struc- 
ture of the density distribution function in a finite grand 
canonical ensemble [7]. Evidently the new approach 
yields estimates of p^(T) and p~{T) of precision ±1-2% 
of Pc or better, for temperatures 1.5 to 2.5 decades closer 
to the critical point. These results confirm convincingly 
that both models belong (as now expected [7,8]) to the 
same (d=3)-dimensional Ising universality class: see be- 
low and the dashed line in Fig. 1. 
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FIG. 1. A log-log plot of the reduced semi-density-jump, 
PS = |[P+(T) - /5"(r)]a3 vs. t={T - Tc)/Tc, where a is 
the hard-sphere diameter, for a HCSW fluid with interac- 
tion range 1.5a (and pc —0.3067) [7] and for the RPM with 
Pc — 0.079 (at a = 5 fine-discretization level [8]). The circles 
report previous estimates for the RPM and HCSW fluid [7] 
employing an equal- weight prescription [10]. The dashed line 
has a slope /3ising = 0.326. 



To outline the established situation [9] , recall that for 
T < Tc the grand canonical equilibrium distribution of 
the density, 7'i(T;p), in a finite system of dimensions 
L'^ with periodic boundary conditions, has two Gaus- 
sian peaks near p^{T) when L^a, where a measures 
the particle size. For large L the two peaks are clearly 
separated and thus provide reasonable estimates for the 
coexisting densities via the equal-weight prescription [10] 
— the open circles in Fig. 1 [7]. However, when Tc is ap- 
proached, finite-size effects, arising from the divergence of 
the correlation length, soon blur the distinction between 
the vapor and liquid states thereby seriously hampering 
the reliable estimation of the coexistence curve. An alter- 
native procedure applicable near Tc is thus imperative. 
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Accordingly, we study the finite-system parameter 
defined [9,11,8] by the dimensionless moment-ratio 

QL{T-{p)L) = {m'')l/{m^)L, m = p-{p)L, (1) 

where denotes a grand canonical expectation value 
at fixed T and fi. As well known, Ql — » | when L — > oo in 
any single-phase region of the (p, T) plane while — > 1 
on the coexistence diameter, pd\an\{T) ^ \ {p^ + At 
criticality, Ql rapidly approaches a universal value Qc 
[8,9,11], e.g., Qc = 0.6236 for d=3 Ising systems. The Q- 
loci, Pq{T] L), on which Ql attains isothermal maxima, 
have recently provided a route to estimating Tc and Pc 
with unprecedented precision [8,12]. 

In the two-phase region it has been known, but little 
appreciated, for some time [10(a), 13], that Ql{T;p) dis- 
plays a surprising singular behavior when L— >oo [14]. 
This is illustrated by the dashed- line plots in Fig. 2, 
which follow directly from the double-peaked structure 
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FIG. 2. Plots of Ql{T;{p)l) vs. p* = {p)a^ for (a) the 
HCSW fluid at T* = 1.200 (< T^* = 1.2182i [12]) and (b) the 
RPM at T* = 0.0500 (<Tc =0.0506.) [8]). The solid lines are 
for (a) L* =L/a = 6,9, 12 and (b) i* =6, 8, 10, 12; the dashed 
lines represent the exact limiting forms for the estimated val- 
ues of p"*" and p~ [15]. 

of Vl{T;p) below [12,13,15] (together with our es- 
timates for /9+ and p~). Specifically, Qoo{T; (p)) ex- 
hibits a discontinuous drop from Qoo ^ f o Qoo = on 
the two-phase boundaries, p~ and and a continuous 
(but nonconvex [12,15]) form for p~<{p)<p'^. For fi- 
nite systems, however, the singularities are rounded and, 
as seen from the histogram-reweighted Monte Carlo sim- 
ulations presented in Fig. 2, Ql{T;{p)l) displays two 
smooth isothermal minima close to p^iT) and p^(T). It 
is notable that while the HCSW data are fairly symmet- 
rical about Pdiam, the RPM displays a remarkably strong 
asymmetry. 

Clearly, it is tempting to extrapolate these minima in 
order to estimate p'^{T) and p~{T) [12]. However, when 



Tc is approached, naive extrapolation fails badly owing 
to the finite-size effects: indeed the graph of Ql{T; (p) l) 
still exhibits two distinct minima at and above T = T^. 
Hence, some more powerful approach is necessary. 

The behavior of Ql{T;{p)l) near criticality can be 
understood via a recently developed 'complete' scaling 
theory that explicitly encompasses pressiire mixing [3 5] . 
Specifically, the full thermodynamics of a one-component 
fluid near the bulk critical point {pc,Tc,Pc) can be de- 
scribed with three relevant scaling fields 

p = p - kot - lofi -\ , 

i=t-hjl- jip-\ , h = p.- kit - j2p-\ , (2) 

where the dimensionless deviations of the pressure 
and chemical potential from criticality are p={p — 
Pc)/pckBTc, and p={p — pc)/k-B.Tc- the coefficients ji 
and j2 measure the degree of pressure mixing, the Yang- 
Yang ratio (« —Tp'^/Cv) being fixed by = —^'2/(1 — 
J2) [3,4]. For a finite box of dimensions L'^ with peri- 
odic boundary conditions, the finite-size scaling hypoth- 
esis now asserts [4,5,12,16] 

p^ppiL-'^Y{x,z), x = DiL^/'', z = Uh/\if, (3) 

where we have used the hyperscaling relation dv = 2 — a 
(valid for d < 4) and, for simplicity, neglected corrections 
to scaling. Note that D and U are nonuniversal ampli- 
tudes (of dimensions L^^/" and L", respectively), while 
Y{x, z) is a universal function that is even in z and in- 
dependent of microscopic details while depending on the 
geometry and the boundary conditions of the system. 
It follows that the full scaling expression for Ql is 

QQ{x,z)[l+AjL-''Qj{x,z)+AiL-^Qi{x,z)+---], (4) 

[12] with exponents and nonuniversal amplitudes 

K = (3/v, Aj=nD^U/pc, 

X={A-l)/v, Ai = {h+n)D^-^/{l-h), (5) 

while the scaling functions Qq , Qj , and Q; depend only 
on derivatives of Y{x,z) thereby being universal. The 
symmetry of Y'(a:, z) implies that Qq is even in z while 
Qj and Qi are odd. Notice that the pressure mix- 
ing coefficient j2 provides the dominant asymmetric in- 
dependent correction (with Ising values K = 0.5l7<A = 
0.895) which, indeed, describes the strong asymmetric 
behavior of Ql{T; {p)l) for the RPM seen in Fig. 2(b). 

Of course, the mean density {p)l also has a scaling 
form which we choose to write as [12] 

y{T-L) = 2[{p)l - Pdiam(T)]/Ap^(T) 

= y[i + AjL-'^yj + AiL-^yi + •••], (6) 

where, again, the scaling functions y{x, z), y, (x, z), and 
yi{x, z) derive from Y{x, z) and are universal, while y is 
odd in z, and yj and 3^; are even. 
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The crucial point here is that pdiam (p"*" + p~) and 
Apoo oc (p"*" — p~ ) embody the desired coexistence val- 
ues p'^{T) and p~(T). Our strategy will be to deter- 
mine values for pdiam and Apoo so that the minima of 
Ql(T; say, Q+{T;L) and Q-{T;L), and their lo- 

cations, p^{T; L) > p^{T; L), satisfy appropriate scaling 
relations. We focus first on Apoc and, to minimize the 
effects of asymmetry (arising from the mixing coefficients 
j2, ji and h), we examine the mean and difference 



5(<9m + <3m), Ay^i„ = i(2/+-y-). (7) 



Now, on evaluating (4) and (6) at 2*;^^ (which asymp- 
totically fixes Qj^j) and formally eliminating x oc tL^/^ be- 
tween the resulting expressions, we see that Qm\n{T;L) 
and A2/niin(r; L) should be related in a way that, to the 
orders displayed, is independent oiT and L and (up to the 
neglected corrections to scaling) reflects only the univer- 
sality class of the critical system under consideration. A 
priori this class is unknown and, indeed, is to be deter- 
mined. However, for any scalar order parameter the two- 
peaked, double-Gaussian structure of Vl{T; p) should be 
reproduced asymptotically when L ^ oo at fixed T <Tc. 
On this basis it is straightforward to calculate the uni- 
versal relation for Qmin^O: we find [12(a)] 



Aymin(«) = 1 



0{q^), g = Q„,i„ln(4/eg„ 



(8) 



which, to this order, is independent of any asymmetry. 

Finally, we can employ our scaling analysis to gen- 
erate the limiting coexistence curve recursively using 
finite-size simulation data for Ql- Appropriate initial 
steps are: (i) Collect data sets {Q^{T;Li),p^{T]Li)} 
for a range of values {ii}" at fixed values of T<Tc. 
(ii) For a value T = To sufficiently low that Qmin^O.03 
[which corresponds to well separated peaks in VlITo] p)], 
choose a density-jump value, say Apxa , independent of 
i, which leads to the best fit of Ay^jj^ = [p+ (Tq; L^) — 
p-{To;Li)]/ApT, vs. q^'^ = q{To; Li) to the relation (8) 
at small q: sec the dashed lines in Fig. 3. In light of 
the scaling relations (4) and (6), the parameter ApT^ 
can then be identified as an estimate for ApadTo). (iii) 
Increase To to Ti=Tq + ATq by a small ATq , chosen 
so that the new set overlaps the previous one. 

(iv) Determine a new value, Apn, so that the plotted 
data display an optimal collapse that extends the previ- 
ous numerical scaling function to larger values of q: see 
the gradual departure of the fits from the dashed lines 
in Fig. 3 as <7 increases. In practice we have found that 
n = 3 distinct box sizes with L3>1.3Li may well suf- 
fice, (v) Repeat steps (iii) and (iv) generating succes- 
sive estimates for Apoo{Tj) for j = 2, 3, ■ • ■. Smaller in- 



crements ATj arc needed as Tj- 



• Tc and the qj^^ increase 



to qc = QJjiin ln(4/e(5Jjji„) (see Fig. 3) so that histogram- 
reweighting procedures are crucial [7,8]. 

Figure 3 presents a scaling plot for the HCSW fiuid 
constructed in this fashion: system sizes L* = L/a = 
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FIG. 3. Scaling plots of (Ay^jn)"^ (for V = 2, 
1/As = 3.O7, and 5) vs. g = Qminln(4/e(5,ni„) for the HCSW 
fluid built up recursively from low q where the dashed lines 
are exact: see Eq. (8). Various symbols, most suppressed for 
clarity, depict results at increasing T,: see text. 



9, 10.5, and 12 were used and led to the estimates shown 
in Fig. 1 for Ap^{T) from |i| ~0.23 down to \t\ ~ 10"*. 
Purely for ease of presentation. Fig. 3 displays (Aymin)"''' 
for selected values of V'- In fact, the scaling analysis indi- 
cates that Ay^niniQ) should diverge like {qc — q)~^ when 
q^qc as T^T^, with qc a imiversal value (depending 
on geometry and boundary conditions) [12(b)]. For the 
HCSW fluid with periodic boundary conditions we find 
QJjjjjj = O.IIO2. To lower precision, the RPM data yield 
the same scaling plots and value of Q^i^ [12(b)]. On the 
other hand, the approximate scaling form proposed by 
Tsypin and Blote [17] for Vl{Tc] p) for ((i=3) Ising mod- 
els gives Q^jjj ~ 0.117, only 6% higher than we observe. 
For {d=2) Ising models we estimate Q^in — 0-28 using 
data in [10]. 

Evidently, the choice of ip = l/P should yield a plot 
that intersects the q axis linearly; indeed, for the Ising 
value, /3is = 0.326, this is so. But, we emphasize that this 
observation plays no role in the calculation of Fig. 1. 

Clearly, uncertainties in choosing ApT^ , Apx^^i , • • • in 
steps (ii) and (iv) will propagate. Well below Tc (where 
care must be taken to ensure two-phase equilibrium) we 
can fit the limiting behavior (8) with a precision of ±1.0% 
or better in Apr/Pc- The overall uncertainties then grow 
by factors of 5-10 as \t\ decreases to 10"'' for the HCSW 
fluid and 10"^ for the RPM: see Fig. 1. 

It is also remarkable that the Apoo{T) estimates im- 
ply values for Tc. For the HCSW fiuid we thus find 
T* ~ 1.21821(2) which lies close to the upper confidence 
limit of the previous estimate T* ~ 1.2179(3) [7]: sec also 
[13(a)] Eq. (5.6). For the RPM we obtain T; c± 0.05069(2) 
which agrees precisely with Ref. [8]. Explicit fits to 
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A/9oc(r), that allow for the leading correction terms, 
yield /3 = 0.324(10) for the HCSW fluid and /? = 0.34(5) 
for the RPM, so providing independent, albeit weaker 
confirmation of the Ising behavior established using data 
confined to T > Tc [7,8]. 

The scaling results (4) and (5) suggest that evidence 
for a pressure-mixing coefficient j2 might be detected in 
finite-size data. Indeed, a detailed calculation [12(b)] of 
the asymmetry seen in the minima of Qt atT = Tc yields 

f^min^^^^) =AjCjL-'^ + AiciL-^ + ..-, (9) 

where Cj and Ci are universal numbers determined by ex- 
pansion coefficients of Y{0, z) about the minima at z^:^^. 
Recall from (5) that Aj is proportional to j2- 

In Fig. 4 we present data for A^j^^{L) for the RPM and 
the HCSW fluid: even by eye, the former strongly suggest 
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FIG. 4. Plots of the critical asymmetry factor A'^i^{L): 
see Eq. (9). The fitted curves use Ising exponent values and 
indicate relatively large pressure mixing in the RPM. 

a leading exponent closer to k = O.5I7 than to A = O.SQe. 
The fits in Fig. 4, using only the two leading terms in (9), 
support this but also indicate a weak j2 contribution of 
opposite sign for the HCSW fluid. Further fairly elabo- 
rate analysis [12] yields j2 = —0.35(7), implying a strong, 
i?^= 0.26(4), Yang- Yang anomaly for the RPM, while 
i2 = 0.042(3) and J?^ = -0.044(3) for the HCSW fluid. 
The latter result is consistent with the earlier, much less 
precise estimate ~ —0.08(12) [7]. 

Finally, to determine the diameter /)diam(r) we com- 
pare ?7min = \{vt, + ym) ^nd Anin(7'; Li). Analysis of the 
two-Gaussian limit [12(b)] yields ymin/^min = \q + 0{q^) 
with q = q — Qmin which is again universal in leading or- 
der. Owing to the asymmetric terms in (4) and (6) the 
analogous scaling plots are now more sensitive to nonuni- 
versal details and exhibit small, L-dependent corrections 



when q approaches qc. Nevertheless, the approach suc- 
ceeds and the critical densities, p*, predicted from the di- 
ameters when Tc are fully consistent with the previ- 
ous, r>Tc estimates [7,8,12]. Details for both the RPM 
and the HCSW fluid will be presented elsewhere [12(b)]. 

In summary, we have shown how the finite-size scal- 
ing information hidden in precise simulation data can be 
systematically extracted via a novel Q-minima recursive 
algorithm to yield coexistence curves far closer to and 
with a much higher precision than previously appeared 
possible. As a byproduct, pressure mixing has been quan- 
titatively resolved. 
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